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Abstract 

A quasi-three-dimensional rotor/statcr analysis has 
been developed for blade-to-blade flows in turbomachin- 
ery. The analysis solves the unsteady Euler or thin- layer 
Navier-Stokes equations in a body- fit ted coordinate sys- 
tem. It accounts for the effects of rotation, radius 
change, and stream- surface thickness. The Baldwin- 
Lomax eddy- viscosity model is used for turbulent flows. 
The equations are integrated in time using a four- stage 
Runge-Kutta scheme with a constant time step. Re- 
sults are shown for the first stage of the Space Shuttle 
Main Engine high pressure fuel turbopump. Euler and 
Navier-Stokes results are compared on the scaled single- 
and multi-passage machine. The method is relatively 
fast and the quasi- three-dimensional formulation is ap- 
plicable to a wide range of turbomachinery geometries. 


Introduction 

The major thrust of the computational analysis of 
turbomachinery to date has been the steady state so- 
lution of isolated blades using mass-averaged inlet and 
exit conditions. Unsteady flows differ from the steady 
solution due to interaction of pressure waves and wakes 
between blade rows. To predict the actual complex flow 
conditions, one must look at the time-accurate solution 
of the entire turbomachine. 

Supersonic computations have been made by M. 
M Rai (1) of unsteady flow in a rotcr/stator config- 
uration. Rai (2) also computed a subsonic flow case 
of a turbine stater using an implicit scheme. He used 
overlaid O- and H- grids around the blade with a con- 
servative treatment at the interface. Gibeling et. al. 
(3) have used an implicit scheme with a distorting grid 
at the interface. Oden and Bass (4) used a finite 
element scheme with an explicit Lax-Wendroff time- 
marching algorithm, along with adaptive mesh refine- 
ment to solve Rai’s supersonic model problem. Giles (5) 
solves the unsteady inviscid Euler equations using an 
explicit Lax-Wfendroff algorithm. A technique is used 
which inclines the computational plane in time such 
that flows can be computed in turbomachines with an 
arbitrary stator/rotor pitch ratio. Three-dimensional 


solutions have been reported by Rai (6). Here the un- 
steady, thin-layer, Navier-Stokes equations are solved 
using a third- order- accurate, implicit, upwind scheme. 
The solution procedure solves for flow in axial- flow tur- 
bomachinery with equal stator/rotcr blade ratios. 

Many of the numerical tools used in the analysis 
of isolated blades can be used fer time-accurate ro- 
tcr/stator interaction analysis. Here the quasi-three- 
dimensional Euler or thin- layer Navier-Stokes equations 
are solved on body fitted C-grids. An explicit four-stage 
Runge-Kutta algorithm is used with a constant time 
step. Computer time is reduced by vectorization. Most 
turbomachinery analyses assume periodicity from blade 
to blade. This makes it possible to analyze one blade 
only. Thrbomachinery is designed with unequal number 
of blades to avoid forced vibration problems. Therefore 
a full unsteady analysis must include at least a few pas- 
sages on each wheel 

The numerical method is presented as well as vis- 
cous results for the first stage of the Space Shuttle Main 
Engine (SSME) high pressure turbopump. Compar- 
isons are made between a scaled 1:1 and 2:3 stator/rotcr 
blade configurations of the turbopump. Inviscid results 
are computed for comparison on the 1:1 stator/rotcr 
configuration. 


Governing Equations 

The axisymmetric (m, 0) coordinate system used 
for the quasi-three-dimensional analysis is shown in 
Fig. 1. Here the m-coordinate is defined by 

+ ( 1 ) 

and the ^-coordinate is defined by: 

0 = (2) 

where O' is fixed in space and 0 rotates with the blade 
row with angular velocity H. The radius r and the 
stream surface thickness h are taken to be known func- 
tions of m. 

The flow equations are written in the (m, 0) system 
(see Refs. 7 and 8) and are then transformed to a gen- 
eral body-fitted (£, rj) system using standard methods. 
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The thin-layer approximation is used to eliminate all 
viscous derivatives in the streamwise (f) direction. The 
final equations are as follows: 


dtq + d(P + a n (G-Re- 1 S)=k (3) 

The in viscid flux terms are given by: 
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The relative contravariant velocity components 
and W* 1 along the £ and rj grid lines are given by: 


' P " 


pV m 

pv 0 r 

e 

P = 7~ 1 


W i = ZmVm + £ 0 W 0 \ W* = rjmVm + T) 0 W 0 (5) 

Here we = ve — rH is the relative tangential velocity 
component. The energy and pressure are given by: 

e = p[C v T+l/2(v 2 m +vl)} (6) 

p=( 7 -l)[ e -l/2(t& + t|)] (7) 

The source term K 2 arises from the centrifugal force 
term in the r-component of the momentum equation. 

K 2 = (pvjj + p - Re~ l <j 22 )r m /r + (p - Re^a^kn/h 

( 8 ) 

where 

, 1 dr 

r ” /r= ;as 

1 dh 
hdm 

Although a similar source term arises from the Coriolis 
force in the ^-momentum equation, the equation has 
been made conservative by multiplying through by r. 
The viscous flux term is 
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where 

5 2 = + *?0 cr 12 

53 = (Vm^ + Vo^y 

S* = (*»m + Vo) d n a 2 + v m S 2 + v 0 S 3 


and the shear stress terms are given by: 

° n = 2nd m v m + XVV 

° 22 = 2p(d 0 v 0 + v m r m )/r + AVF 

0-33 = 2 piv m hm/h + AV ■ V 

(712 = n{d m V6 - v 0 r m /r + l/rd g v m ) 

2 (10) 
A W = -- nldmVm + v m (r m /r + hm/h) + l/rd 0 v 0 ] 

Using the thin-layer approximation, the shear stress 
terms are evaluated by replacing d m with rj m d n and 
(l/r)d 0 with fj 0 drj. Here a = y/j p/p is the sonic vebc- 
ity, and the normalized thermal conductivity k equals 
one. 

The equations are nondimens ionalized by arbitrary 
reference quantities (here the inlet total density and 
critical sonic velocity define the reference state), and 
the Reynolds number Re and the Prandtl number Pr 
must be specified in terms of that state. These equa- 
tions assume that the specific heats C p and C v and the 
Prandtl number are constant, that Stokes’ hypothesis 
A = —2/3/2 is valid, and that the effective viscosity may 
be written as 


M — fMamdnar IMurbulent 

The transformation metrics are found using 

(t %Y j {X ~%) < u > 

where the Jacobian is given by 

J = {m i 6 n -m n 0 ( )- 1 (12) 

Overbars in Eqs. (49) denote a rescaling of the metrics: 

le = Ze/ r > Ve = W r i T 1 = rhJ- 1 (13) 

For turbulent flows the two-layer eddy-viscosity 
model developed by Baldwin and Lomax (9) is used. 
In the (m, 6) coordinate system the wall shear r w and 
vorticity w required by the model are given by 

r w = <?i 2 w = Mi^mVe + 1 /rd 0 v m - verm/r) w 

w = l(dmVo - 1 /rd 0 v m + v 0 r m /r) 

(14) 


Computational Grid 
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Body-fitted C-type grids for this work were gen- 
erated using the GRAPE code developed by Sorenson 


(10,11). Figure 2 shows typical grids around the first 
stage stater and rote** of the SSME. The stator grid 
was modified by adding a single grid line parallel to the 
exit to allow approximately one cell of overlap with the 
rotor grid. The rotor grid was also selectively refined 
by doubling the number of £=constant grid lines in the 
inlet region. Both grids were modified by adding one 
line of dummy points (not shown) corresponding to in- 
terior points from neighboring passages. These points 
are used for imposing the periodic or overlap boundary 
conditions. 


Multistage Runge-Kutta Algorithm 


An explicit multistage Runge-Kutta algorithm 
based on the work of Jameson (12,13) is used to ad- 
vance the flow equations in time. Given the residual R 
from a central-difference representation of Eq. (3), a 
/c-stage scheme may be written as follows: 

q(°)=q n 

q W = q(o) - ai JA t\Bq {0) ~ D(q^)} 


?<*> = q ( o) _ ak JAt\Rq( k -V -D(g( 0 >)] 

g n+1 = q( k ) 

(15) 

The time step is given by At = min( At,-,y) where 
A tij is calculated from an inviscid stability analysis of 
Eqs. (3-15) to be 

A* . . = CFL 

+ M le + ay/l?n + le + 

m W 

where l m and 1$ are reciprocal length scales given by 
= | | + | tfm | 

(17) 

In this work we have used a four-stage scheme with the 
standard coefficients Oj=( 1/4, 1/3, 1/2, 1). This scheme 
is second- order accurate in time and has a Courant limit 
CFL= 2.8. The results shown here were calculated with 
a maximum Courant number of 2.5, which occurs in 
fine-grid regions near the wall. Courant numbers in the 
coarse-grid cere flow regions are typically two orders of 
magnitude smaller for viscous flows. 


Artificial Dissipation 


Dissipative terms consisting of fourth and second 
differences are added to prevent odd-even point decou- 
pling and to allow shock capturing respectively. The 
dissipative term D in Eq. (15) is given by: 

Dq = (D i +D n )q (18) 

The ^-direction operator is given by: 

D^q = C(V 2 q(t - (19) 


where 

»^A tij 

is a coefficient that partially cancells similar terms in 
Eq. (15). We have found that using the spatially- 
varying time step A Uj (Eq. (16)) as a coefficient in 
the dissipation is much less dissipative than using the 
minimum time step. 

The terms V 2 and V 4 are given by: 


V 2 =M 2 max(i/ i+ 1 j, u id , n-i,y) 
V 4 = max(0, fj .4 — V2) 


where 


and 


P»-t-i,y ~ + Pt— i,y 

P<+i,j+2pi,y + Pi-i, j 


( 20 ) 

( 21 ) 


M2 = 0(1) 

P4 = 

( 22 ) 

In smooth regions of the flow the dissipative terms 
are of third order and thus do not detract from the for- 
mal second- order accuracy of the scheme. Near shocks 
i/ij is large and the second-difference dissipation be- 
comes locally of first order. 


Boundary Conditions 


Inlet 


At the inlet to the stator, total pressure, total 
temperature T, and the whirl rvo are specified. At 
each time step the upstream-running Riemann invari- 
ant R~ = v m — 2a/ (7 — 1) is extrapolated to the inlet. 
The axial velocity component v m is found using 
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V m - 


Interface Formulation 



Density and energy are found using isentropic relations. 
This is a non-reflective approximation to the axial mo- 
mentum equation that is first- order accurate in space 
but zero- order in time. However in the present re- 
sults we see virtually no unsteady behavior at the inlet 
boundary. 

Outlet 

At the exit from the rotor, static pressure is spec- 
ified and the other flow quantities are found using 
second-order extrapolation of p, pv m) and pv$. The 
interaction between the specified exit pressure and the 
inlet conditions sets the mass flow through the machine. 
This exit condition is strictly reflective and may influ- 
ence the unsteady solution. We hope to sort out the 
effects of this condition in future work. 

Solid Wall 

Blade surface pressures are found from the normal 
momentum equation: 

(tmVm + + {riL + f}j)d n p 

= -pW^^md^Vm-^fjgd^ve) +pv e (r] m V 0 -7j e v m )r m /r 

(24) 

Surface velocities are found from the tangency or no-slip 
conditions for inviscid or viscous flows, and the blade 
temperature is specified. 

Periodic Boundaries 

The periodic boundaries are solved in the same way 
as the interior points. The grid fcr each blade is over- 
laped one grid cell at its periodic boundary such that 
the grid points are coincident with the grid points of an 
adjacent blade. These overlap dummy grid points are 
updated when an adjacent blade is solved and are then 
used in the next time step integration. For equal pitch 
blade rows the flow variables on the spatially periodic 
boundaries are equal on each blade. This is not true 
for multi-passage calculations. Here, periodic boundary 
conditions are enforced only on boundaries of blades 
that are encompassed by the blade ratio. Flow vari- 
ables on interior boundaries are updated using quanti- 
ties from adjacent blade solutions. 


For a given displacement between stator and rotcr 
blade rows, the solution on the interface must be deter- 
mined. Figure 3 shows an enlargement of the interface 
geometry of the stater and rotor grids. The upstream 
boundary of the rotor is coincident with the £ — 2 and 
f = M — 1 grid lines in the stater grid. The stator grid 
extends into the rotor grid such that the exit boundary 
grid points of the stator lie within the first cell of the 
rotcr domain. 

The solution is updated at the interface by inter- 
polating the flux variables in the stator computational 
domain to obtain rotor quantities and interpolating the 
flux variables in the rotor computational domain to ob- 
tain stator quantities. The interface is updated in this 
way after every Runge-Kutta integration sweep of the 
stator and rotor domains. Rai (1,2,6) has stressed the 
importance of conservative treatment of the interface. 
Our present formulation is nonconservative and the im- 
portance is unknown. 


Solution Procedure 
Initial Conditions 

Given a desired leading edge velocity triangle and 
exit flow angle fcr each blade, the initial conditions 
are calculated using an analytic solution of the one- 
dimensional flow equations (see Ref. 7). These ini- 
tial conditions do not especially speed the convergence 
to a periodic unsteady solution, but they do provide a 
smooth transition between the inlet and exit boundary 
conditions. 

Startup 

Since the initial conditions are not known, it is de- 
sirable to reach the conditions where the initial tran- 
sdents are no longer in the solution domain. From the 
initial analytic flow condi tions, the solution is integrated 
using a local time step. This allows the mass flow to 
stabilize and the wake region to develop quickly. A con- 
stant minimum time step is then used to compute the 
time accurate solution. 

Data Management 

In a steady state analysis, the solution can be con- 
verged rapidly by using a local time step and multi-grid. 


T 


4 



When a time- accurate solution is desired, a multigrid 
procedure cannot be used and a minimum time step 
must be used in the integration procedure. This is fur- 
ther restricted when rotor /stator interaction is consid- 
ered, since the solution must be integrated using a con- 
stant minimum time step based on the computational 
domain of two blades. r 

Our current strategy is to update each blade row a 
full time step (four Runge-Kutta stages), then to update 
the inflow, exit, blade, and interface boundaries. The 
solution is considered converged when the loading on 
the blades becomes periodic in time. 

The solution for each passage is stored as a sepa- 
rate array in core and accessed as the solver alternates 
from one blade to the next. The flow solver is rela- 
tively unaffected by the number of passages so that the 
multi- passage problem becomes a problem of data man- 
agement. 


Results 

The solution procedure has been applied to the first 
stage turbine rotor of the SSME fuel fcurbopump. The 
actual SSME stator /rotor blade count ratio for the first 
stage is 41:63. The numerical procedure was first ap- 
plied to a single passage blade row configuration where 
the stator was scaled to the rotor pitch. The scaling 
was done such that the pitch-to-chord ratio remained 
unchanged. 

The grids generated using the GRAPE code have 
viscous s pacings of 0.0002 in. away from the blade. The 
stator grid has 115 x 3 1 points. The rotor grid has 197 x 
41 grid points. It was necessary to cluster grid points 
at the inlet of the rotor to get the necessary resolution 
to capture the wake from the upstream stator. 

The flow Mach number varies from . 148 at the star 
ter inlet through about .472 between the blade rows to 
.201 at the rotor exit. The inlet absolute flow angle to 
the stator is 0.0 degrees and the relative inlet flow angle 
to the rotor is -31.0 degrees. Using the above grid spac- 
ing, flow conditions and a maximum Corn* ant number 
of 2.5, the time step is calculated based on the inviscid 
stability limit for the four-stage Runge-Kutta scheme. 
The time step is such that it takes 11262 steps to move 
one rotor pitch. This takes about 18 minutes on the 
Cray X-MP for the equal pitch computation. A solu- 
tion is considered converged when the unsteady lift on 
each blade beccxnes periodic. For the equal passage mar 
chine this convergence criterion is met after seven rotor 
blade passings. This requires 2.10 hours of computer 
time. 


The unsteady lift diagrams fer the single passage 
machine are shown in Fig. 4. Fluctuations in lift on the 
stator are minimal but are more dramatic on the rotor. 
The minimum lift takes place when the rotor begins to 
pass through the wake of the stator. 

The static pressure envelope for the stator, Fig. 
5, indicates that the only influence of the downstream 
rotor is on the uncovered portion of the suction surface 
of the stator. The major portion of the blade remains 
unchanged from the time-averaged solution. The static 
pressure envelope of the rotor, Fig. 6, shows a larger 
influence of the stater wake at the leading edge of the 
rotor blade which influences the rest of the blade. The 
adverse pressure gradient seen on the pressure surface 
of the blade forces the large cove separation seen in Fig. 
7. 

Figures 8 through 11 show Mach contours as the 
rotor moves through the stator wakes. Remnants of the 
wakes from other stator blades can be seen in the rota- 
passage. When entropy contours are plotted in Fig. 12, 
this is more apparent. 

Next we present results for a multi-passage configu- 
ration with a 2:3 blade count that is close to that of the 
actual machine. The unsteady lift diagrams are shown 
in Fig. 13 for a single stafca and a single rotor. The 
affect of the larger upstream stators on the downstream 
rotas is much more dramatic than seen in the single 
passage computation. Again the stator row sees only 
small fluctuations due to the rotor passing by down- 
stream. The three fluctuations seen in the lift show the 
influence of the three passing downstream rota's over a 
periodic pitch. A rotor passes through two stator wakes 
over one periodic pitch rotation. The rota sees large 
fluctuations in its lift, with the minimum lift occurring 
when the rota’s leading edge encounters the wake of the 
upstream stator. A periodic solution is obtained after 
about seven rota* blade passings. This takes about six 
hours of computer time and requires about 1.25 times 
the storage of a single-passage computation. 

The static pressure envelope for the multi-passage 
stator, Fig. 14, shows the fluctuation in the pressure 
only on the uncovered pation of the suctiai surface of 
the blade. The pressure gradient is less adverse than 
that seen on the single passage machine. The static 
pressure envelope on the rota, Fig. 15, shows a larger 
static pressure fluctuation along the entire surface of the 
blade. The adverse pressure gradient seen on the pres- 
sure surface near the leading edge generates a smaller 
cove separation on the multi-passage machine, Fig. 16. 

Figures 17 through 20 show Mach contours fa the 
multi- passage machine. The asymmetry of the flow is 
apparent in the rotor passages. Here the stator wakes 
can be seen at different locations of rotation. Again the 
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entropy contours, Fig 21, clearly show the asymmetry 
of the different rotor passages as the wakes move down- 
stream. 

The invisdd solution of a single- passage machine 
was computed by matching the static pressure at the 
exit of the rotor with that of the viscous solution. The 
grids used in this calculation allowed a larger time step 
to be used such that 1874 iterations were necessary to 
move one pitch rotation. 

The static pressure envelope on the stator, Fig. 22, 
shows fluctuations similar to those in the viscous com- 
putations. On the rotor, Fig. 23, the static pressure 
plots show a larger adverse pressure gradient at the lead- 
ing edge and a faster recovery. Also, the fluctuations in 
pressure along the rotor blade due to the upstream sta- 
ter in the viscous solution are not as prominent in the 
inviscid solution. 


Concluding Remarks 

The quasi-t hree-dimensional Euler or the thin- 
layer, Navier-Sfcokes equations are solved for unsteady 
turbomachinery flows. These equations are written in 
general coordinates for an axisymmetric stream surface, 
and account for the effects of blade row rotation, radius 
change and stream- surface thickness. 

A four-stage Runge-Kutta scheme based on the 
work of Jameson is used to predict time-accurate re- 
sults. A nonconservative interface formulation and 
other data management techniques allow the solution 
of rotor /stater interaction problems in both single and 
multi- passage machines. 

Viscous results on the first stage of the Space Shut- 
tle Main Engine high pressure fuel turbopump have 
shown the analysis to be viable. Results have been 
shown on the scaled single- as well as the multi-passage 
machine. The fact that the analysis is quasi-t lire e- 
dimensional will allow the time-accurate solution of ra r 
dial flow turbomachinery. 
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Figure 1. - Quasi-three-dimensional stream surface and coordinate 
system for a centrifugal compressor. 



Figure 2. - Computational grids about SSME. 



Figure 3. - Enlargement of grid at the interface. 
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Figure 4. - Unsteady lift diagram for 1:1. 


Figure 5. - Static pressure envelope on stator (1:1 viscous) 
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Figure 6. - Static pressure envelope on rotor (1:1 viscous). 


Figure 7. - Relative velocity vectors on rotor showing 
cove separation (1:1). 
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Figure 8. - Mach number contours (1:1, j pitch rotation). 


Figure 9. - Mach number contours (1:1, | pitch rotation). 




Figure 10. • Mach number contours (1:1, | pitch rotation). 


Figure 11. • Mach number contours (1:1, aligned). 
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Figure 17. - Mach number contours (2:3, aligned). 


Figure 18. - Mach number contours (2:3, J rotor pitch 
rotation). 



Figure 19. - Mach number contours (2:3, | rotor pitch 
rotation). 


Figure 20. - Mach number contours (2:3, | rotor pitch 
rotation). 
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